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ABSTRACT 


This  study  determines  how  to  invest  limited  resources  to  increase  the  resilience  of  an 
infrastructure  system  against  both  non-deliberate  and  deliberate  events.  We  propose  an 
optimization  model  that  seeks  the  best  defensive  investment  for  a  weighted  combination 
of  deliberate  and  non-deliberate  events.  We  formulate  the  general  problem  and  conduct 
numerical  analysis  using  a  specific  infrastructure  system  as  a  concrete  example.  We 
perform  parametric  analysis  on  the  combined  model  in  order  to  explore  the  way  in  which 
different  solutions  depend  on  the  distributed  weights  and  yield  insight  into  the  best 
investment  decisions. 
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EXECUTIVE  SUMMARY 


Owners  and  operators  of  critical  infrastructure  systems  face  a  challenge  when  it  comes  to 
investing  limited  resources  to  maintain  continuity  of  function  in  the  presence  of 
disruptive  events.  On  the  one  hand,  they  need  to  protect  against  the  most  commonly 
occurring  disruptions  (e.g.,  failures  or  weather  events).  On  the  other  hand,  they  need  to 
protect  against  worst-case  disruptions,  as  might  be  caused  by  an  intelligent  adversary 
who  uses  insider  knowledge  to  deliberately  interrupt  system  function.  Often  these  two 
types  of  potential  disruptions  point  to  different  parts  of  the  system  as  “critical,”  and  in 
that  case  this  tension  can  lead  to  a  dilemma  in  terms  of  what  to  protect. 

In  this  study,  we  consider  how  to  invest  limited  resources  to  increase  the 
resilience  of  an  infrastructure  system  against  both  non-deliberate  and  deliberate  events. 
We  represent  the  function  of  the  infrastructure  as  a  network  flow  problem  and  use  an 
optimization  model  that  seeks  the  best  defensive  investment  for  a  weighted  combination 
of  deliberate  and  non-deliberate  events.  We  conduct  parametric  analysis  on  the  combined 
model  in  order  to  explore  the  way  in  which  different  solutions  depend  on  the  distributed 
weights  and  yield  insight  into  the  best  investment  decisions.  From  our  analysis,  we 
discover  that  it  is  possible  to  choose  the  best  defense  for  that  specific  network  by 
mapping  the  frontier  between  the  worst-case  and  random  disruptions.  The  results  of  that 
frontier  gives  us  insight  into  (1)  the  minimum  number  of  edges  to  protect,  (2)  the 
maximum  defense  for  both  worst-case  and  random  disruptions,  and  (3)  a  moderate  cost- 
effective  defense  for  both  worst-case  and  random  disruptions 

Our  investigation  of  a  small  network  system  provides  proof-of-concept  that  the 
technique  we  propose  can  be  effective,  but  there  is  more  work  to  be  done.  In  particular, 
our  solutions  are  computationally  intensive  even  for  small  systems,  and  additional 
research  to  reduce  the  required  computational  effort  is  needed  before  the  technique  can 
be  used  at  large  scale. 
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I.  INTRODUCTION 


The  United  States  (U.S.)  has  critical  infrastructure  systems  (e.g.,  electric  power, 
transportation,  telecommunications,  water)  that  need  to  be  protected  from  disruptive 
events  that  are  non-deliberate  (e.g.,  natural  disasters,  accidents,  failures)  and  deliberate 
(e.g.,  vandalism,  sabotage,  terrorism).  Recently,  these  systems  have  been  the  targets 
terrorist  attacks  (September  11,  2001)  and  natural  disasters  (Hurricane  Katrina  in  2005 
and  Hurricane  Sandy  in  2012)  which  caused  billions  of  dollars  in  damage.  In  response  to 
these  and  other  disruptive  events,  the  U.S.  is  spending  billions  of  tax  payer  dollars 
annually  (e.g.,  $50  billion  in  2014)  to  strengthen  infrastructure  (The  White  House  2014). 
Of  particular  interest  is  how  to  allocate  these  funds  in  a  manner  that  increases  the 
resilience  of  our  critical  infrastructure  systems. 

In  the  reliability  engineering  literature  (e.g.,  Hwang  et  al.  1981  and  Billinton  and 
Allan  1992)  and  the  risk  literature  (e.g.,  Bedford  and  Cook  2001),  non-deliberate 
disruptions  are  commonly  modeled  as  random  events  using  probabilities.  However, 
deliberate  disruptions  are  more  appropriately  represented  as  worst-case  events 
using  game  theory  (e.g.,  Brown  et  al.  2006).  These  two  perspectives  often  point  to 
different  parts  of  an  infrastructure  system  as  “critical.”  This  leaves  policy  makers  in  a 
quandary  about  how  to  invest  a  limited  defensive  budget  to  make  the  infrastructure 
system  more  resilient. 

The  goal  of  this  study  is  to  determine  how  to  invest  limited  resources  to  increase 
the  resilience  of  an  infrastructure  against  both  non-deliberate  and  deliberate  events. 
This  infonnation  will  provide  decision  makers  with  insight  into  which  assets  to  protect 
and  how  much  emphasis  to  put  toward  preparing  for  natural  disasters  vice  a 
terrorist  attack. 

We  propose  an  optimization  model  that  seeks  the  best  defensive  investment 
for  a  weighted  combination  of  deliberate  and  non-deliberate  events.  We  formulate  the 
general  problem  and  conduct  numerical  analysis  using  a  specific  infrastructure  system 
as  a  concrete  example.  We  perform  parametric  analysis  on  the  weights  in  the  combined 
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model  in  order  to  explore  the  way  in  which  different  solutions  depend  on  the  distributed 
weights  and  yield  insight  into  the  best  investment  decisions. 
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II.  LITERATURE  REVIEW 


Presidential  Decision  Directive  21  (PPD21),  issued  in  February  2013,  establishes 
national  policy  on  critical  infrastructure  security  and  resilience  (The  White  House  2013). 
Following  the  definition  established  previously  in  the  USA  PATRIOT  Act  (Department 
of  Justice  (DOJ)  2001),  PPD21  defines  critical  infrastructure  to  include  “systems  and 
assets,  whether  physical  or  virtual,  so  vital  to  the  United  States  that  the  incapacity  or 
destruction  of  such  systems  and  assets  would  have  a  debilitating  impact  on  security, 
national  economic  security,  national  public  health  or  safety,  or  any  combination  of  those 
matters.” 

PPD21  specifically  states  “It  is  the  policy  of  the  United  States  to  strengthen  the 
security  and  resilience  of  its  critical  infrastructure...  to  take  proactive  steps  to  manage 
risk  and  strengthen  the  security  and  resilience  of  the  Nation's  critical  infrastructure, 
considering  all  hazards  that  could  have  a  debilitating  impact  on  national  security, 
economic  stability,  public  health  and  safety,  or  any  combination  thereof.  These  efforts 
shall  seek  to  reduce  vulnerabilities,  minimize  consequences,  identify  and  disrupt  threats, 
and  hasten  response  and  recovery  efforts  related  to  critical  infrastructure.”  In  its 
definition  of  “all  hazards,”  PPD2 1  refers  to  “a  threat  or  an  incident,  natural  or  manmade 
[that] . . .  includes  natural  disasters,  cyber  incidents,  industrial  accidents,  pandemics,  acts 
of  terrorism,  sabotage,  and  destructive  criminal  activity  targeting  critical  infrastructure.” 

In  response  to  PPD21,  the  Department  of  Homeland  Security  (DHS)  issued  an 
updated  National  Infrastructure  Protection  Plan  (NIPP)  in  December  2013.  The  2013 
edition  of  the  NIPP  “continues  to  focus  on  risk  management  as  the  foundation  of  critical 
infrastructure  security  and  resilience”  (p.  4). 

The  application  of  risk  analysis  to  critical  infrastructure  systems  follows  a  long 
history  in  the  application  of  probabilistic  risk  assessment  (PRA)  to  complex  engineering 
systems  (see  Bedford  and  Cook  2001  for  an  introduction,  or  National  Aeronautics  and 
Space  Administration  (NASA)  2011  for  a  recent  application  of  PRA  to  space  mission 
planning).  Following  the  attacks  of  9/11,  much  of  the  work  on  risk  to  critical 
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infrastructure  has  focused  on  the  application  of  PRA  to  terrorism  risk  (e.g.,  Pate-Cornell 
and  Guikema  2002,  Garrick  et  al.  2004,  Parnell,  et  al.  2005,  Willis  2007).  The 
Department  of  Homeland  Security  has  promoted  the  use  of  probabilistic  techniques  when 
assessing  the  risk  critical  infrastructure  systems  (DHS  2013). 

A  different  line  of  research  has  focused  on  the  ability  of  an  intelligent  adversary 
to  deliberately  disrupt  system  function.  Specifically,  there  is  now  considerable  work  on 
the  use  of  system  interdiction  models  for  understanding  homeland  security  problems  and 
the  risk  to  critical  infrastructure  systems — these  models  are  known  as  Attacker-Defender 
(AD)  and  Defender-Attacker-Defender  (DAD)  problems  (Brown  et  al.,  2005,  Brown  et 
al.,  2006,  Alderson  et  al.,  2014a,  b).  Applications  of  this  model  have  been  used  in  a 
variety  of  systems:  fuel  (Ileto  2011),  cargo  (De  La  Cruz  2011),  the  electric  grid 
(Sahneron  et  al.  2004,  2009),  and  military  networks  (Burton  2013,  Long  2013). 

In  this  thesis,  we  propose  a  model  to  defend  our  infrastructure  against  both  non- 
deliberate  hazards  and  deliberate  threats.  We  use  a  linear  combination  of  the  two 
techniques  previously  described  in  order  to  map  the  efficient  frontier  of  the  tradeoffs 
between  the  consequences  of  a  worst-case  attack  and  expected  consequence  over  a  set  of 
random  events  that  will  give  decision  makers  the  ability  to  choose  the  type  of  protection 
needed  for  his  or  her  particular  network. 

Previous  work  by  Bier  (2007),  Zhuang  and  Bier  (2007),  and  Hausken  et  al.  (2009) 
has  looked  at  the  challenge  of  choosing  what  to  protect,  particularly  when  facing  “all 
hazards”  to  include  deliberate  and  non-deliberate  events.  However,  these  studies  have 
abstracted  the  consequence  of  the  event  to  a  simple  known  function.  As  demonstrated  by 
Alderson  et  al.  (2013),  in  general  it  is  difficult  to  approximate  the  consequences  of  losing 
one  or  more  infrastructure  components,  or  to  pick  the  most  vital  components  in  a  network 
system.  An  explicit  objective  of  this  research  is  to  connect  the  operational  models  of 
infrastructure  function  and  consequence  assessment  with  the  need  to  identify  protective 
measures  that  work  for  both  deliberate  and  non-deliberate  disruptive  events. 
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III.  MODEL  DEVELOPMENT 


In  this  chapter,  we  present  a  sequence  of  models  that  we  use  to  perform  our 
analysis.  Throughout,  we  follow  the  concepts  and  notation  in  Alderson  et  al.  (2014a).  A 
specific  implementation  of  these  models  follows  in  Chapter  IV. 

A.  NORMAL  OPERATION 

In  this  section,  we  provide  an  overview  of  the  sequence  of  models  to  follow.  Our 
starting  point  is  an  Operator  Model  of  the  form 

min  /  (w,x,y)  ,  (1) 

_pe7(w) 

where  w  denotes  the  fixed  design  of  the  system,  x  denotes  the  operational  setting  (e.g., 
the  status  of  each  component  of  the  system),  and  Y(w)  denotes  the  set  of  feasible 
activities  available  to  the  system  operator,  given  design  w. 

The  Operator  Model  represents  the  decision  making  of  an  entity,  called  the 
operator,  who  chooses  the  activities  in  the  system.  The  Operator  Model  is  prescriptive;  its 
solution  y*  =argmin  f(w,x,y)  represents  the  activities  in  the  system  that  yield  the 

yey(w) 

minimum  operating  cost. 

The  Operator  Model  (1)  is  a  mathematical  program  that  can  be  solved  using 
various  techniques  (see  Alderson  et  al.  2014a  for  a  detailed  discussion). 

B.  WORST-CASE  DISRUPTIONS 

We  are  interested  in  understanding  how  well  the  system  performs  in  the  presence 
of  worst-case  disruptions  and/or  most  likely  disruptions.  For  a  fixed  design  w,  we  define 

P(w)  =  max  min  / (w,x,y) ,  (2) 

xgX  7g7(w) 

where  X  is  the  set  of  feasible  disruptive  events,  and  x  =argmax  min  f(w,x,y )  is  a 

xeX  yzY( w) 

worst-case  disruptive  event.  Here,  P(w)  represents  the  performance  of  the  system  in  the 

presence  of  the  worst-case  disruptive  event.  We  refer  to  (2)  as  the  Attacker  Model;  it 

builds  on  the  previous  fonnulation  but  has  additional  elements. 
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In  general,  P(w)  can  be  evaluated  three  ways: 

1 .  by  exhaustive  enumeration  of  x  eX , 

2.  by  taking  the  dual  of  the  inner  Operator  Model  to  transform  the  max-min 
of  the  Attacker  Model  into  a  max-max  problem  (e.g.,  Brown  et  ah,  2006), 
and 

3.  using  Bender’s  decomposition  (e.g.,  Benders,  1962). 

Given  a  set  of  possible  designs  W ,  our  goal  is  to  find  a  design  weW  that  yields 
“good”  perfonnance  in  the  presence  of  disruptive  events. 

If  we  are  only  concerned  with  defending  against  worst-case  disruptions,  then  we 
solve  a  problem  of  the  form:  min  P(w) .  This  is  simply  a  detenninistic  Defender- 

weW 

Attacker-Defender  (DAD)  problem.  Let  W*D  a  W  denote  the  set  of  optimal  designs  that 
minimize  P(w).  That  is,  w*  e  W*D  <=>  w*  =  arg  min  P(w)  . 

weW 

In  principle,  solving  the  Defender  Model  requires  nothing  more  than  enumerating 
every  possible  combination  of  defense  and  attack,  then  solving  the  corresponding 
Operator  Model  for  each,  and  then  finding  the  one  that  yields  the  lowest  cost.  This  is  the 
first  half  of  our  combined  model  which  will  determine  the  effect  of  the  worst  case.  The 
second  half  of  the  model  is  the  probabilistic  model  that  optimizes  defense  against 
hazards. 

C.  PROBABILISTIC  DISRUPTION  (MOST  LIKELY) 

For  a  fixed  design  w,  define 

Q(w)  =E-  min  f{w,x,y),  (3) 

y&Y  (w) 

where  x  is  a  random  variable  having  a  specified  distribution,  and  E.  denotes  the 

expectation  relative  to  x .  Here,  Q(w)  represents  the  average  performance  of  the  system 
in  the  presence  of  a  random  disruptive  event. 

In  general,  Q(w)  can  be  evaluated  two  ways: 

1 .  via  Monte  Carlo  sampling  on  x ,  and 
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2.  using  exhaustive  enumeration  if  x  is  discrete  to  compute  an  exact 
expectation. 


If  we  are  only  concerned  with  defending  against  random  disruptive  events,  then 
we  solve  a  problem  of  the  form:  min  Q(w) .  This  is  simply  a  stochastic  optimization 

weW 

problem.  Let  W*  eW  denote  the  set  of  optimal  designs  that  minimize  Q(w).  That  is, 
w*  e  W*  o  w*  =  argmin  Q(w)  . 

weW 


D.  DEFENDING  AGAINST  WORST-CASE  AND  MOST  LIKELY 

In  practice,  we  are  concerned  with  both  random  disruptive  events  and  worst-case 
disruptive  events.  To  determine  a  good  design,  we  must  take  into  considerations  the 
following  two  possibilities. 

1.  Observe  that  if  WD  nWs*  ^  0  ,  then  we  have  an  optimal  solution  to  both 
problems  and  there  is  no  issue.  We  simply  choose  w*  e  (W'r)  n  IT  * )  . 

2.  However,  in  general  we  expect  W*D  =  0  .  More  strongly,  our 
experience  is  often  that  the  worst-case  disruptive  event  is  very  different 
from  the  average  case  disruptive  event.  This  creates  a  dilemma  in  terms  of 
where  to  invest  defensive  resources  to  make  the  system  more  resilient. 

We  propose  to  study  the  combined  problem  as 

Z  *  (Z)  =  min[/L P(w)  +  (1  -  Z)(2(w)]  where  0  <  A  <  1 .  (4) 

weW 

Note  that  when  A  —  0 ,  Eq.(4)  corresponds  to  the  stochastic  average-case.  In 
contrast,  when  A  —  1 ,  Eq.(4)  corresponds  to  the  deterministic  worst-case. 


7 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


8 


IV.  CASE  STUDY 


In  this  thesis,  we  use  a  small  numerical  example  to  illustrate  our  technique. 
Following  Alderson  et  al.  (2014b),  Figure  1  displays  a  notional  (fuel  distribution) 
infrastructure  system  made  of  nodes  and  links.  The  system  is  a  single  commodity 
transshipment  model,  designed  to  move  fuel  from  supply  locations  to  demand  locations. 


Figure  1.  Notional  fuel  distribution  network  (from  Alderson  et  al.  2014b). 


In  Figure  1,  the  white  circles  represent  locations  with  demand  =  1  and  the  black 
circles  represent  locations  with  supply  =10.  All  links  are  bidirectional  and  have  per  unit 
flow  usage  cost  =  1.  Each  link  has  capacity  =15  and  the  penalty  for  unsatisfied  demand 
per  node  =10.  Usage  penalty  for  an  interdicted  arc  =  20  per  unit  of  flow  (e.g.,  hiring  a 
truck  to  drive  the  unit  of  flow  to  bypass  the  link).  The  objective  is  to  minimize  the  sum  of 
the  flow  cost  for  fuel  deliveries  and  any  penalties  from  unsatisfied  demands. 

As  shown  in  Figure  1,  Nodes  3,  4,  and  16  have  two  (parallel,  redundant) 
connections  to  the  rest  of  the  network.  This  network  has  been  built  to  be  N- 1  reliable, 
meaning  that  the  loss  of  any  single  link  does  not  disconnect  any  node. 


9 


A.  OPERATOR  MODEL 


We  formulate  the  Operator  Model  for  this  system  as  follows. 

Indices  and  Sets 

n  EN  nodes  (alias  i,j) 

[i,j\  EE  undirected  edge  between  nodes  i  and  j,  i<j 

(/,/)  E  A  directed  arc  from  node  i  to  node  j 

[ij]  EE  (/</)  A  (( ij )  E A  A  (j,i )  E A)) 

dED  defenses,  £>={0,1} :  0  represents  undefended,  1  represents 

defended  (and  invulnerable) 

Data  [units] 

Cy  per  unit  cost  of  traversing  arc  (ij)  EA  [dollars/barrel] 

Ujj  upper  bound  on  total  (undirected)  flow  on  edge  [i,j]  EE 

[barrels] 


Xij  1  if  edge  [i,j]  EE  damaged,  0  otherwise  [binary] 

qc'.  per  unit  penalty  cost  of  traversing  arc  (i,j)  EA  if  damaged 
[dollars/barrel] 

bn  fuel  supply  at  node  n  EN  [barrels] 

(-demand  for  bn  <0  ) 

vn  per  unit  penalty  cost  for  demand  shortfall  n  EN 

[dollars/barrel] 

w''  1  if  defense  option  d  is  implemented  on  edge  [i,j]  EE, 

0  otherwise  [binary] 
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Decision  Variables  [units] 


Yy  flow  on  arc  (ij)  E  A  for  defense  d  [barrels] 

Sn  fuel  shortfall  at  node  n  EN  [barrels] 


Formulation: 


min 

Y,S 

X  X  [k + v'l-% )  YS + (c> + vpv )  Yi 

deD  [i,j]<=E 

]+XvA 

ueN 

(DO) 

s.t. 

y  y  Yd~y  y  Yd-s  <t 

deD  j:(n,j)eA  dEDi:(i,n)EA 

Vn  e  N 

k]  (CD 

Yd  +Yj.  <  U..Wy 

IJ  Jl  ij  IJ 

V[/,  j]  <=A,deD 

[k]  (°2) 

o 

Al 

V[z,  j]  eA,d  <eD 

(D3) 

Vn  eN 

(D4) 

Discussion 

The  objective  function  (DO)  combines  the  total  flow  cost  and  the  total  penalty 
cost.  Constraints  (Dl)  enforce  balance  of  flow  at  each  node.  Stipulations  (D2)  and  (D3) 
ensure  bounds  on  decision  variables.  This  formulation  implements  cost-based 
interdiction — that  is,  damage  to  an  arc  makes  it  extremely  expensive  but  not  infeasible — 
which  makes  the  problem  easier  to  solve  computationally.  In  the  above  example,  we  have 
bs=bio=  10  and  bn=- 1  otherwise.  In  addition,  we  set  cy=l ,  q\  =  0,  V(i,  j)  e  A 

Uy  =  100,  V(z,y)  e  A  ,q®  =  10,V(z',y)  e  A ,  Uij=  15  for  all  [ij]  EE,  and  vn  =10  for  all  nE N. 
There  are  no  arcs  defended  in  the  base  case;  therefore,  =0  and  w'i  =1  for  all  [ij]  EE. 

Solving  this  problem,  we  obtain  a  minimum-cost  solution  consisting  of  flows  that 
meet  all  the  demands  (see  Figure  2),  with  the  flows  labeled  in  blue.  Our  findings  are 
consistent  to  the  base  solution  findings  in  Alderson  et  al.  (2014b).  We  are  concerned  with 
the  ability  of  the  system  to  perfonn  in  the  presence  of  either  random  (probabilistic)  or 
deliberate  (worst-case)  disruptive  events. 
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Figure  2.  Operator  Model  base  solution 


B.  ATTACKER  MODEL 

For  this  notional  fuel  system,  we  model  the  attacker’s  problem  as  follows. 


Additional  Data[units] 

ry  “cost”  to  break  edge  [/,/]  E  E  [cardinality] 

attack_budget  budget  constraint  on  the  number  of  simultaneous  attacks 

[cardinality] 


Additional  Decision  Variables  [units] 

Xy  1  if  attacker  breaks  edge  [/,/]  EE,  0  otherwise  [binary] 


Formulation 
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(ADO) 


P(w)  =  max  min  i  z  [(W*x+(<v+«;w]+ Z”A 

JeZ)  [ij]e£  weiV 

s.t.  (Dl),  (D2),  (D3)(D4) 

y  njXij  <  attack  _  budget 

UJ]^E 

A.ejO,!}  V[f,y]  e  E 


(ADI) 

(AD2) 


The  objective  function  (ADO)  is  the  same  as  that  for  the  Operator  Model  (DO), 
except  that  parameters  Xy  have  been  replaced  by  decision  variables  Xy.  Constraint  (ADI) 
limits  the  number  of  simultaneous  attacks,  and  the  cost  to  attack  each  arc  can  be  different. 
Stimpulations  (AD2)  require  that  attacks  are  binary.  We  note  that  qy  =  0  implies  that  arc 
(ij)  is  effectively  invulnerable,  because  attacking  it  does  not  increase  the  flow  cost  for  the 
operator. 

In  the  above  example,  we  model  parallel  arcs  as  costing  twice  as  much  to  attack. 
That  is,  we  have  r2j  =  r4$=  r  12,16  =  2  and  all  other  ry  =  1. 

The  solution  for  the  Attacker  Model  depends  on  the  attacker  budget  (the  number 
of  attacks  that  the  attacker  can  afford).  The  attacker  seeks  to  maximize  the  cost  to  the 
operator.  In  Figure  3,  an  attack  on  edge  [10,13]  represents  the  worst  case  attack  to  the 
operator  when  the  attacker  can  afford  to  target  only  a  single  edge.  The  operator  reroutes 
his  flow  as  shown  in  blue  and  the  overall  operator  cost  increases  from  25  to  33.  An 
attacker  having  the  ability  to  target  a  large  number  of  edges  could  render  our  network 
practically  useless  as  shown  in  Figure  4.  Here,  the  attacker  has  the  budget  to  target  five 
edges  and  focuses  on  the  arcs  around  the  supply  nodes  to  cripple  flow  to  the  entire 
network  with  exception  to  node  4.  The  attacker  picks  the  best  place  to  attack  depending 
on  how  many  attacks  allotted  by  the  attacking  budget.  These  results  are  also  consistent 
with  Alderson  et  al.  (2014b). 
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Figure  3. 


Figure  4. 


Operator  Model  solution  with  5  attacks 


C.  DEFENDING  AGAINST  A  WORST-CASE  ATTACK 

In  this  thesis,  we  restrict  defensive  actions  to  the  “hardening”  of  edges,  such  that 
targeting  an  edge  does  not  increase  its  usage  cost.  We  formulate  the  following  model. 

Additional  Data[units] 

defense _budget  budget  constraint  on  the  number  of  defenses  [cardinality] 
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Additional  Decision  Variables  [units] 


w: 


1  if  defense  option  d  is  implemented  on  edge  [ij]  EE, 
0  otherwise  [binary] 


Formulation 


min  max  min  X  £  [(^  +  <^)^+(c.i 

'  dcD  [i,j]eE 

+4W?]+I>A 

nsN 

(DADO) 

s.t.  (D1)(AD1)(AD2)(D3)(D4) 

Yij  +  Y*  <  UqW* 

V[/,  j]  eE,d  eD 

(DAD\) 

^  W(J  <  defense _budget 

[Ej]&E 

(DAD2) 

II 

t'xb 

V[z, ./]  e  E 

(DAD3) 

TT  o-i} 

V[i,  /]  e  E,d  eD 

(DAD4) 

The  objective  (DADO)  includes  the  cost  of  flow  over  existing  (and  possibly 
damaged  or  protected)  edges  and  penalties  for  unmet  demand.  Contraints  (DAD1) 
ensures  that  flow  on  edge  [i,j]  does  not  exceed  flow  capacity  Uy.  The  constraint  (DAD2) 
requires  that  the  cost  of  all  defenses  fall  within  the  existing  defense  budget  using  the 
binary  variable  W? .  Constraints  (DAD3)  ensure  that  flow  on  each  edge  \ij\  is  on  either 

the  protected  or  unprotected  edge  and  not  both.  Here,  the  pair  IT1  =  1  and  IT'1  =0  denotes  a 

defended  edge  and  vice  versa  for  undefended.  Stipulations  (DAD4)  require  binary 
defense  decisions. 

In  order  to  be  consistent  with  Alderson  et  al.  (2014b)  we  will  give  examples  of 
defense  with  an  attacker  budget  of  three  attacks  as  shown  in  Figure  5.  For  the  purposes  of 
this  study  we  only  focus  on  hardening  edges. 
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Figure  5.  A.  AD  Model:  worst-case  attack  involving  three  edges.  In  the  absence 
of  any  defenses,  edges  [2,7],  [10,13],  and  [11,15]  are  targeted.  B.  DAD 
Model:  the  corresponding  worst-case  attack  with  one  edge  protected. 

Edge  [10,13]  is  defended,  and  edges  [6,10],  [7,8],  and  [9,13]  are 
targeted. 

The  defender  model  protects  edge  [10,  13],  note  that  it  is  the  same  arc  the  attacker 
chose  as  the  one  arc  attack  in  Figure  3.  The  attacker  picks  another  set  of  arcs  to  attack  in 
order  to  maximize  cost  for  the  operator  that  does  not  include  the  defended  arc.  With  edge 
[10,13]  defended  and  in  the  presence  of  three  targeted  edges,  the  operator  is  able  to 
satisfy  demand  at  seven  nodes,  instead  of  only  six  with  no  arcs  defended. 

The  defender  can  limit  the  attacker  by  strengthening  more  arcs  at  shown  in  Figure 
6.  However,  we  cannot  defend  all  arcs  because  of  budget  constraints.  Therefore, 
decisions  about  what  arcs  to  protect  must  be  made  within  those  constraints  to  minimize 
cost  to  the  operator. 
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Figure  6.  Worst-case  disruption  with  3  edges  targeted  and  5  edges  protected. 

D.  EXPECTATION  MODEL 

We  use  the  symbol  x  to  represent  a  random  event  in  place  of  x ,  which  denotes  a 
set  of  edges  that  have  been  damaged.  We  calculate  the  expected  cost  due  to  random 
disruptions,  and  will  not  focus  on  a  single,  worst-case  event.  For  example, 

H  =  H2>,  <  2 1  represents  the  set  of  all  events  with  two  or  fewer  links  broken.  For  any 

element  xeS,  let  P~  be  the  probability  that  x  occurs,  i.e.,  { 0 < P~  <  1  Vx e 3  , 

V  R  =1  }.  For  simplicity,  we  assume  the  probability  that  an  edge  will  go  down  is 

xeH 

independent  between  edges,  giving  it  a  binomial  distribution  with  x  becoming  a  random 
variable.  We  can  think  of  E  as  a  list  of  the  enumerated  (all)  possible  attacks  with 
associated  probabilities. 

Additional  Sets 

£  e  E  set  of  (enumerated)  outcomes 

Additional  data  [units] 
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x-c  c:  {0, 1}^  set  of  binary  vectors,  one  element  for  each  edge,  where  Jc f.  =  1 

if  edge  [z',y']e£'is  damaged,  and  =0  if  the  edge  remains  intact,  in 
outcome 


Formulation 

Q(w)  =  E4(F(w,^))  ,  where 

F(w,c)  =  minX  Z  [(%  +  qfi! W?  +  (cji  +  44 )*?  ]  +  Z  VA 

flfeD[(,y]6£  «e./V 

si.  (D3)(D4) 

Z  Z  ^-Z  Z  Y'-SnZb*  VneN 

d&D  j:(n,j)eA  d&D  i:(i,n)eA 

Y?  +  Yjf  <Uijwfj  V[z,y]e  A,d  eh 


(E0) 


[4]  (El) 
[>f  ]  (E2) 


Whereas  the  AD  model  yields  both  the  worst-case  attack  and  associated  best- 
response  in  the  objective  value,  the  result  of  the  ED  model  is  solely  the  objective  value. 

Defending  against  the  average  case  disruption  means  we  must  take  into  account 
the  arcs  most  likely  to  fail  and  either  strengthen  or  build  around  those  arcs  in  each  x  . 

E.  DEFENDING  AGAINST  THE  EXPECTED  CASE 

Just  like  in  the  DAD  model,  we  are  defending  edges  in  order  to  minimize  the 
system  cost  in  the  presence  of  disruption  by  carefully  choosing  which  edges  to  protect. 

The  possible  defenses  will  be  the  same  as  in  the  DAD  model,  just  like  in  Figure  7. 
Flowever,  the  result  will  be  the  objective  value  comprised  of  the  products  from  this 
design  and  all  possible  failures  in  E  creating  the  defender-expectation-defender  (DED) 
case. 
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Figure  7.  Example  Defense  in  which  edges  [2,7],  [7,8],  [9,13],  [10,1 1],  and 
[10,13]  are  defended. 


Formulation 

In  order  to  determine  the  best  defense  for  the  expected  value  problem,  we  use 
Benders  decomposition  (Benders,  1962),  an  iterative  algorithm  in  which  we  use  the  index 
k  to  denote  an  iteration,  and  K  to  denote  the  current  iteration  of  the  algorithm.  The 
subproblem  for  any  fixed  defense  w  and  any  specific  outcome  ^  e  H  is  just  F  ( w,  , 

and  can  be  solved  as  a  linear  program  with  associated  duals  )  as  given  in  (El) 

and  (E2).  Each  of  these  dual  solutions  provides  a  linear  lower  bound  on  the  objective 
function  value  of  F(»,£)  as  a  function  of  the  defense: 

^(•.£E2>X+  I 

neN  [i,j]&E  d<=D 

The  weighted  sum  of  these  linear  functions  is  thus  a  lower  bound  on  the  expected 
cost  as  a  function  of  defense,  so  we  then  define: 

K=YjP^n 

and 

to  be  the  weighted  combinations  of  the  optimal  dual  solutions  to  each  of  the  subproblems 
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F  [wk  solved  in  iteration  k,  yielding  a  single  cut  for  the  set  of  constraints  (SMI)  in 
the  master  problem.  We  therefore  have  the  following  stochastic  master  problem: 


min  Z«„ 

(SMO) 

s.t. 

neN  [i,j]<EE  deD 

k  =  \,...,K 

(SMI) 

IX =1 

deD 

V[/,y]eE 

(SM2) 

y,  Wy  <  defense  budget 

u,n*E 

(SM3) 

The  objective  (SMO)  represents  the  expected  cost  of  operating  in  the  face  of 
random  events.  We  minimize  this  objective  through  each  iteration  k  as  shown  in  (SMI). 
Each  edge  [i,j]  has  two  possibilities,  denoted  with  d,  and  only  one  can  be  used  at  a  time 
shown  in  constraint  (SM2).  The  constraint  (SM3)  ensures  that  the  defenses  used  are 
limited  to  defense _budget.  We  describe  the  procedure  for  solving  this  problem  in  the 
context  of  the  combined  model,  which  now  follows. 

F.  COMBINED  MODEL 

If  we  define  vvfloch  to  be  the  weight  given  to  the  stochastic  (i.e.,  expected- value) 
objective,  where  0  <  wtstoch  <  1 ,  and  wtdet  =  1  -  vvflocli  to  be  the  remaining  weight  placed  on 

the  deterministic  (i.e.,  worst-case)  objective,  then  we  can  formulate  the  following 
combined  master  problem: 
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min 

z,w 


(CMO) 


s.t. 


^^stoch^stoch  ^^deC^det 


(SM1),(SM2),(SM3) 


Z<le,  ^  Z  Z  +(C7/+^y)y/V_ 


rfeZ)  [ij]e£ 


-Z  Vnsn 


k  =  l,...,K 


neN 


(CM2) 

(CM3) 


Z  Z  y*-Z  Z  ri?-S‘.<b„VneN,k=l,...,K 

deD  j:[n,j]eA  deDi-{i,n]eA 

\/[i,j\eE,dtD,k  =  \,...,K  (CM4) 
V (i,j)  £  A,d  e  D,k  =  (CM5) 

Sf>0  Vn  e  N,k  =  (CM6) 


Y  +Y  <u~W‘ 
y  p  u  v 


Y*k>  0 

~ik 


Here  the  defense  decisions,  W,  are  made  before  any  realization  of  an  attack  or  a  random 
failure,  and  the  evaluation  of  the  defense  chosen  at  iteration  A:  is  a  convex  combination  of 
the  evaluation  of  each  of  the  two  separate  subproblems. 

At  each  iteration  K  of  Benders  decomposition,  our  algorithm  solves  each  of  the 
two  subproblems  (stochastic  and  deterministic)  to  optimality  for  a  fixed  defense,  wK . 
The  constraints  from  the  stochastic  master  problem  transfer  directly  to  this  combined 
master  problem.  The  detenninistic  subproblem  yields  a  single  worst-case  attack,  xK  ,  and 
this  provides  a  lower  bound  on  the  cost  of  the  optimal  response  for  that  particular  attack, 
although  only  in  terms  of  the  flow  variables  in  response  to  that  attack,  and  this 
contributes  a  new  cut  to  the  set  of  constraints  (CM2)  in  the  master  problem.  We  therefore 
must  also  create  a  block  of  nonnegative  variables,  YK ,  to  model  that  response,  and  two 
additional  blocks  of  constraints  to  ensure  they  (a)  conform  to  balance  of  flow  (CM3),  and 
thus  act  like  admissible  flows  in  our  network,  and  (b)  correspond  appropriately  to  any 
defense,  W  (CM4).  The  algorithm  then  solves  the  master  problem  for  iteration  K,  and 
generates  a  new  defense,  wk  1 .  After  the  subproblems  have  been  solved,  we  have  an 
accurate  evaluation  of  the  combined  cost  of  the  current  defense  plan,  \vK  .  If  this  cost  is 
lower  than  any  seen  so  far,  we  record  wK  as  a  new  incumbent  solution  and  take  its  cost 
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as  the  new  upper  bound.  Each  solve  of  the  master  problem  provides  a  lower  bound  on  the 
optimal  solution,  and  these  values  are  monotonically  nondecreasing.  If,  at  any  point,  the 
relative  difference  between  the  best  upper  bound  derived  from  the  subproblems  and  the 
best  lower  bound  from  the  master  problem  is  less  than  a  predetermined  optimality  gap,  s, 
the  algorithm  terminates  and  reports  the  current  incumbent  as  an  ^-optimal  solution.  This 
combined  model  uses  the  duals  from  the  stochastic  model  and  combines  them  with  the 
worst  case  model  to  produce  cuts. 
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V.  RESULTS  AND  ANALYSIS 


In  this  chapter,  we  explore  the  differences  between  the  “worst-case”  and  “most 
likely”  disruptions  for  the  notional  fuel  distribution  network  in  Figure  1.  We  use  the 
combined  model  to  perfonn  our  analysis,  where  in  equation  (CMO)  we  have  convex 
weights  wtdet=  a  and  wtstoch  =  1  -  a. 


A.  ANALYSIS  OF  FIVE  SCENARIOS 

We  ran  the  combined  model  for  five  scenarios,  each  described  in  the  subsections 

below. 


1.  Two  Attacks  and  One  Defense 

Figure  8  displays  the  results  of  the  combined  model  for  our  notional  fuel  system 
in  the  case  of  two  attacks  and  one  defense.  For  different  values  of  a,  we  solve  the 
combined  model  to  within  an  optimality  tolerance  of  0.001.  On  a  computer  with  a 
dual-core  1.2  GHz  processor,  the  solution  time  for  all  scenarios  involving  two 
attacks  was  approximately  2  hours  for  all  eleven  cases  (a=0.0  through  a=1.0).  For  each 
value  of  a,  Figure  8  shows  the  value  for  Zstoch  and  Zdet  along  with  their  combination 
weighted  by  a . 

When  defending  solely  against  the  worst-case  attack  (a=l),  the  best  defense  is 
[2,7]  and  the  worst  two-link  attack  is  [7,8]  and  [8,12].  When  defending  solely  against 
random  events  (a=0),  the  best  defense  is  also  [2,7].  Note  that  the  value  Zdet  is 
significantly  higher  than  Zstoch,  meaning  that  the  worst-case  disruption  is  significantly 
worse  than  the  average  case.  Because  the  worst-case  solution  does  not  vary  with  a,  the 
values  for  Zstof/iand  Zdet  are  insensitive  to  changes  in  a.  As  a  result,  the  combined 
objective  value  is  not  meaningful  on  its  own,  rather  it  represents  simply  the  weighted 
combination  of  these  two  constant  values. 
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Weight  (a) 


-Z_det 

Zstoch 

Objective_(CMO) 


Figure  8.  The  objective  values  for  two  attacks  and  one  defense  versus  the  weight 
for  given  weight  a. 


2.  Two  Attacks  and  Two  Defenses 

Just  as  in  the  previous  case,  the  worst-case  solution  does  not  vary  with  a,  the 
values  for  Zstoch  and  Zdet  are  insensitive  to  changes  in  a.  When  defending  solely  against 
the  worst-case  attack  (a=l),  the  best  defense  is  [2,7]  and  [7,8];  while,  the  worst  two-link 
attack  is  [1,2]  and  [9,13].  When  defending  solely  against  random  events  (a=0),  the  best 
defense  is  also  [2,7]  and  [7,8]. 

Table  1  summarizes  the  results  from  the  scenarios  involving  no  more  than  two 
broken  links. 


Table  1.  Optimal  links  to  defend  for  different  defensive  budgets  for  cases 
when  there  are  no  more  than  two  broken  links.  For  this  small 
network,  we  observe  that  the  single  best  link  to  protect  is  [2,7],  and 
when  defending  two  links  the  best  two  are  [2,7]  and  [7,8]. 


Stoch 

Combined 

Det 

Def 

(a  =0) 

(a  =  5) 

(a=l) 

1 

[2,7] 

[2,7] 

[2,7] 

2 

[2,7] [7,8] 

r2,7]r7,8] 

[2,7][7,8] 
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3. 


Three  Attacks  and  One  Defense 


The  results  of  the  combined  model  for  our  notional  fuel  system  in  the  case  of 
three  attacks  and  one  defense  has  a  shape  similar  to  the  two-attack  case.  For  this  scenario, 
we  solve  the  combined  model  for  different  values  of  a  to  within  an  optimality  tolerance 
of  0.1.  On  a  computer  with  a  dual-core  1.2  GHz  processor,  the  solution  time  for  was 
approximately  6  hours  for  all  eleven  cases  (a=0.0  through  a=1.0).  Table  2  displays  the 
results  of  this  computation. 


Table  2.  Combined  model  results  for  three  attacks  and  one  defense 


a 

edge  to  defend 

worst  attack 

objective  (CM0) 

obj  Zdet 

obj  Zstoch 

0.0  (Stoch) 

12,71 

19,13116,10117,81 

29.47955282 

80 

29.47955282 

0.1 

T2,71 

r9,13H6, 10H7, 81 

34.53159754 

80 

29.47955282 

0.2 

T2,71 

[9,131[6,10U7,81 

39.58364226 

80 

29.47955282 

0.3 

T2,71 

[9,131[6,101[7,81 

44.63568697 

80 

29.47955282 

0.4 

T2,71 

[9,13116,10117,81 

49.68773169 

80 

29.47955282 

0.5 

[2,7] 

[9,13][6,10][7,8] 

54.73977641 

80 

29.47955282 

0.6 

T2,71 

[9,13116,10117,81 

59.79182113 

80 

29.47955282 

0.7 

12,71 

[9,131[6,101[7,81 

64.84386585 

80 

29.47955282 

0.8 

T2,71 

[9,131[6,101[7,81 

69.89591056 

80 

29.47955282 

0.9 

n3,ioi 

[9,131[6,101[7,81 

75.02040504 

80 

30.2040504 

1.0  (Del) 

[13,101 

[9,13][6,10][7,8] 

80 

80 

30.2040504 

When  defending  solely  against  the  worst-case  attack  (a=l),  the  best  defense  is 
[13,10]  and  the  worst  three-link  attack  is  [7,8],  [9,13]  and  [6,10]  as  shown  in  Table  2. 
While  defending  solely  against  random  events  (a=0),  the  best  defense  is  [2,7].  There  is  a 
difference  in  the  solution  for  the  worst-case  and  the  solution  for  random  events;  however, 
the  solution  doesn’t  change  until  a  >  0.9  toward  the  worst-case.  The  difference  in  cost  to 
the  operator  between  protecting  edge  [2,7]  and  [13,10]  is  one.  Similar  to  the  two  attack 
scenarios,  the  value  Zdet  is  significantly  higher  than  Zstoch.  There  is  more  value  in 

protecting  [2,7]  for  the  stochastic  case  than  [13,10],  but  the  values  for  Zsloch  are  similar  for 
the  two  cases.  Overall,  the  values  for  Zstoch  and  Zdet  are  insensitive  to  changes  in  a,  but 
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are  far  from  each  other.  Again,  the  combined  objective  value  represents  simply  the 
weighted  combination  of  these  two  constant  values. 


4.  Three  Attacks  and  Two  Defenses 

We  similarly  solve  for  the  results  of  the  combined  model  for  our  notional  fuel 
system  in  the  case  of  three  attacks  and  two  defenses.  On  a  computer  with  a  dual-core  1.2 
GHz  processor,  the  solution  time  for  was  approximately  10  hours  for  all  eleven  cases 
(a=0.0  through  a=1.0).  Table  3  displays  the  results  of  this  computation. 


Table  3.  Combined  model  results  for  three  attacks  and  two  defenses 


a 

edge  to  defend 

worst  attack 

objective  (CM0) 

obj  Zdet 

obj  Z stoch 

0.0  (Stoch) 

12,7117,81 

[1,21[13, 101111, 15] 

28.03713818 

72 

28.03713818 

0.1 

[2,7][7,8] 

[1,2][13, 10][11, 15] 

32.43342436 

72 

28.03713818 

0.2 

113,10117,81 

[13, 14][8, 12][10, 11] 

36.47357923 

67 

28.84197404 

0.3 

113,10117,81 

[13,141[8,121[10,111 

40.28938182 

67 

28.84197404 

0.4 

[13,10][7,8] 

[13, 14][8, 12][10, 11] 

44.10518442 

67 

28.84197404 

0.5 

[13,10][7,8] 

[13, 14][8, 12][10, 11] 

47.92098702 

67 

28.84197404 

0.6 

[13,10117,8] 

[13,141[8,121[10,111 

51.73678961 

67 

28.84197404 

0.7 

[13,10][7,8] 

[13, 14][8, 12][10, 11] 

55.55259221 

67 

28.84197404 

0.8 

[13,10][7,8] 

[13, 14][8, 12][10, 11] 

59.36839481 

67 

28.84197404 

0.9 

ri3, 10117, 81 

[13,141[8,121[10,111 

63.1841974 

67 

28.84197404 

1.0  (Del) 

113,10117,81 

[13, 14][8, 12][10, 11] 

67 

67 

28.84197404 

The  difference  between  defending  solely  against  random  events  (a=0)  and  the 
worst-case  attack  (a=l)  is  the  defending  of  edge  [2,7]  vice  [13,10].  Although  [13,10]  and 
[2,7]  yield  the  two  highest  operating  costs  for  attacking  a  single  edge,  it  is  more 
important  to  protect  only  one  of  them  along  with  [7,8]  due  the  construction  of  the 
network  as  shown  in  Figure  1.  Unlike  the  previous  defense  for  three  attacks,  the  solution 
turns  to  the  worst-case  when  a  >0.2,  meaning,  when  coupled  with  the  defense  of  [7,8]  it 
is  more  beneficial  for  both  models  to  protect  [13,10]  vice  [2,7]. 

5.  Three  Attacks  and  Three  Defenses 

We  similarly  solve  for  the  results  of  the  combined  model  for  our  notional  fuel 
system  in  the  case  of  three  attacks  and  three  defenses.  On  a  computer  with  a  dual-core  1 .2 
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GHz  processor,  the  solution  time  for  was  approximately  13  hours  for  all  eleven  cases 
(a=0.0  through  a=1.0).  Table  4  displays  the  results  of  this  computation. 


Table  4.  Combined  model  results  for  three  attacks  and  three  defenses 


a 

edge  to  defend 

worst  attack 

objective  (CM0) 

obj  Zdet 

obj  Zstoch 

0.0  (Stoch) 

r2, 7H7, 8H9, 131 

[1,21[13,101[11,151 

27.45686501 

72 

27.45686501 

0.1 

12, 7H7, 8U13, 101 

[8,121[13,141[10,111 

31.26976353 

67 

27.29973726 

0.2 

[2,7][7,8][13,10] 

[8, 12][13, 14][10, 11] 

35.23978981 

67 

27.29973726 

0.3 

113,10117,8118,121 

[2,7119,131110,111 

39.12518556 

64 

28.4645508 

0.4 

113,10117,8118,121 

[2,71[9,131[10,111 

42.67873048 

64 

28.4645508 

0.5 

[13,101[7,8][8,12] 

[2,7][9,13][10,11] 

46.2322754 

64 

28.4645508 

0.6 

[13,101[7,81[8,121 

[2, 71[9, 13U10, 111 

49.78582032 

64 

28.4645508 

0.7 

[13,10][7,8][8,12] 

[2,7][9,13][10,11] 

53.33936524 

64 

28.4645508 

0.8 

[13,101[7,81[8,121 

R, 7H9, 13H10, 111 

56.89291016 

64 

28.4645508 

0.9 

[13,101[7,81[8,121 

[2,71[9,131[10,111 

60.44645508 

64 

28.4645508 

1.0  (Det) 

[13,10][7,8][10,1 1] 

[2,7][9,13][8,12] 

64 

64 

28.68990798 

The  three  attack  and  three  defense  solution  for  the  combined  model  isn’t  as 
intuitive  as  the  previous  two  scenarios.  When  0.3  <  a  <  0.9  ,  the  optimal  solution  defends 
edges  [8,12],  [13,10],  and  [7,8],  Showing  that  moderate  protection  against  both  random 
and  worst-case  disruptions  is  separate  from  the  solutions  for  either  event  as  shown  in 
Table  2. 

Figure  9  shows  that  as  the  values  for  the  separate  models  change,  the  objective 
function  still  maintains  the  linear  characteristic  as  expected;  however  the  operator  costs 
for  Zdetgoes  down  for  each  change  in  defense  while  Zstoch  goes  up  slightly  for  the 

increase  in  a.  The  worst  case  and  expectation  models  both  chose  the  same  edges  to 
defend  for  the  two  attack  scenario  as  shown  in  the  Table  1.  This  result  makes  it  easy  for  a 
decision  maker  to  choose  which  items  to  defend  for  the  best  possible  protection.  The 
problem  for  decision  makers  is  shown  in  the  case  of  three  attacks. 
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■Z_det 
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Objective_(CMO) 


Figure  9.  The  objective  values  for  three  attacks  and  three  defenses  versus  the 
weight  for  given  weight  a. 


B.  DISCUSSION 

In  the  case  of  random  disruptions,  with  a  limited  number  of  defenses,  we  might 
not  be  able  to  protect  our  system  against  the  (potentially  enormous)  set  of  disruptions 
near  the  “center”  of  the  distribution.  However,  if  there  are  only  a  few  disruptions  that 
cause  significant  damage  (i.e.,  well  above  the  mean)  then  we  can  apply  our  limited 
resources  to  reduce  their  impact  and  therefore  have  a  noticeable  effect  on  the  expected 
cost  over  the  entire  distribution  of  random  disruptions. 

Figure  10  displays  the  overall  value  of  having  one,  two,  or  three  defenses,  when 
protecting  against  three  attacks.  There  are  diminishing  returns  for  the  third  defense  used. 
To  save  money  the  decision  maker  should  use  a  two  defense  approach  vice  three  in  this 
case. 
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Weight  (a) 


'0 


defense1 


defense1 

defense1 


'3 


defense1 


Figure  10.  Multiple  defense  options  against  the  three  attacks  scenario:  Objective 
value  for  given  defense  with  weight  a  in  equation  (CMO),  where  wtdet 
=  a  and  wtst0ch  =  1  -  a. 
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VI.  CONCLUSIONS 


In  this  thesis,  we  have  presented  a  model  for  defensive  investment  that  combines 
the  threat  of  a  worst-case  disruption  with  the  risk  of  a  random  disruption.  Our 
combination  model  has  the  potential  to  help  decision  makers  come  up  with  a  good 
defense  that  fits  their  particular  network.  From  our  analysis  we  discover  that  it  is  possible 
to  choose  the  best  defense  for  that  given  network  by  mapping  the  frontier  between  the 
worst-case  and  random  disruptions. 

Our  investigation  of  a  small  network  system  provides  proof-of-concept  that  the 
technique  we  propose  can  be  effective,  but  there  is  more  work  to  be  done.  We  discovered 
that  the  stochastic  model  runs  longer  than  the  combined  model.  Interestingly,  this 
suggests  that  knowing  how  to  defend  against  the  worst-case  potentially  provides  insight 
into  how  to  defend  against  random  events.  However,  even  for  the  small  example  here,  the 
program  takes  a  long  period  of  time  to  solve  the  combined  model  for  three  attacks  and 
two  or  more  defenses.  Therefore,  additional  research  to  reduce  the  required 
computational  effort  is  needed  before  the  technique  can  be  used  at  large  scale. 
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APPENDIX 


Our  computational  models  and  solution  were  obtained  with  the  Python 
programming  language  and  the  Cooper-Pyomo  optimization  module  using  CPLEX  as  the 
solver.  The  program  is  solved  using  five  modules.  The  combo_main .  py  module 
contains  the  outer  loop  for  multiple  runs  of  the  model  and  is  the  driver  that  branches  into 
the  other  modules.  We  have  a  separate  module,  phat_code .  py,  to  create  the 
distribution  for  the  probability  of  damaging  one  or  more  arcs.  The  network  is  built  using 
the  networkx  module,  in  netx_graph  .  py,  by  building  the  list  of  edges  and  nodes.  The 
file  combo_model .  py  creates  the  pyomo  models  for  the  master  and  sub  problems  used 
to  create  cuts  and  implement  benders.  The  file  pyomo_model  .py  creates  and  solves 
the  minimum  cost  flow  model.  The  code  for  the  modules  are  listed  below: 


#Begin  combo_main .py 

import  netx_graph 

import  pyo_model 

import  combo_model 

import  phat_code 

import  matplotlib . pyplot  as  pit 

#open  file 

f ile=  open ( ' ComboResults . txt '  , '  w '  ) 

file . write (' Weight  (Worst),  Arcs  Defended,  Objective  value.  Worst  attack, Worst 
value, Random  value\n') 

test  =  [ 0 , 0 . 1 , 0 . 2 , 0 . 3 , 0 . 4 , 0 . 5 , 0 . 6 , 0 . 7 , 0 . 8 , 0 . 9 , 1 ] 
for  weight  in  test: 

if  _ name _  ==  ' _ main _ '  : 

#establish  nodes,  edges,  start_nodes,  and  double_edges 

#This  will  eventually  occur  in  a  subroutine  that  parses  input  files 

nodes  =  ['I',  '2',  '3',  '4',  '5',  '6',  '7',  '9',  'll',  '12',  '13',  '14',  '15',  '16'] 

edges= 

[('1',  '2'),  ('2',  '7'),  ('1',  '5'),  ('5',  '9'),  ('9',  '13'),  ('13',  '10'),  ('13',  '14'),  (’1 
4',  '15'),  ('6',  '10'),  ('6',  '7'),  ('7',  '8'),  ('8',  '12'),  ('10',  'll'),  ('11',  '12'),  ('ll 
'15')  ] 


start_nodes  =  ('10', '8') 

double_edges  =[('2','3'),('8','4'),('12','16')] 

#build  graph 

Grph,  vuln,  double_nodes  =  netx_graph . newmodel (nodes,  start_nodes, 
edges,  double_edges) 
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#establish  single-edge  probabilities  and  initialize  xhat 

phat,  xhat  =  phat_code . phat_edges (Grph,  vuln) 

defenses  =  ['O',  '  1 '  ] 

w_hat  =  { } 

for  (i,j)  in  vuln: 

w_hat [i, j , ' 0 ' ]  =  1 
w_hat [i, j , ' 1 ' ]  =  0 


#Create  dictionaries  for  node  supplies  and  vulnerable  edge  capacities 
supply  =  { } 

for  item  in  Grph . nodes () : 
if  item  in  start_nodes: 

supply[item]  =  10 
elif  item  in  double_nodes : 

supply [item]  =  0 
else : 

supply[item]  =  -1 
capacity  =  { } 
for  (i,j)  in  vuln: 

capacity [ i, j ]  =  10 

#generate  probability  map  for  all  events  with  two  or  fewer  failures 
P  =  phat_code . create_Prob (Grph, vuln, phat, 3) 
print  'Number  of  events:  '  +  str(len(P)) 
max_iterations  =  50 

nodepi  =  {}  #dictionary  of  node  duals,  by  iteration 
edgemu  =  {}  #dictionary  of  arc  capacity  duals,  by  iteration 
xhats  =  { }  #dictionary  of  worst-case  attacks,  by  iteration 


LB=0  #  can  set  this  based  on  first  MP  solve... 

UB=10000  #  should  come  up  with  a  better  way  of  handling  this 

master_tol  =0.10  #  relative  gap  tolerance  for  terminating  benders 

det_wt  =  weight 
stoch_wt  =  1-weight 

for  iteration  in  range (max_iterations ) : 
expected_cost  =  0 
cumul_p  =  0 

for  ( scennum, x_tilde)  in  enumerate (P . keys ()) : 

xhat  =  phat_code . generate_xhat_from_x_tilde (vuln, x_tilde, xhat ) 
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model,  instance 

pyo_model .build_model (Grph, xhat, w_hat, vuln, defenses, supply, capacity) 

results  =  pyo_model . solve_model ( instance) 

instance . load (results) 

obj ective_value  =  instance . FlowCost ( ) 

if  scennum  ==  0: 

worst_scen  =  x_tilde 

worst_obj  =  obj ective_value 

for  node  in  Grph . nodes () : 

dualval 

instance . dual . getValue ( instance . FlowBalance [node] ) 

if  abs (dualval)  >  le-08: 

nodepi [node, iteration]  =  P [x_tilde] *dualval 
else : 

nodepi [node, iteration]  =  0.0 
for  (i,j)  in  vuln: 

for  defense  in  defenses: 
dualval= 

instance . dual . getValue ( instance . EdgeCapacity [ i, j , defense] ) 

if  abs (dualval)  >  le-08: 

edgemu [i, j , defense, iteration] = 


P [x_tilde ] * dualval 


else : 

edgemu [i, j , defense, iteration]  =  0.0 


else : 

if  worst_obj  <  ob j ective_value : 
worst_obj  =  obj ective_value 
worst_scen  =  x_tilde 
for  node  in  Grph . nodes () : 
dualval= 

instance. dual. getValue (instance . FlowBalance [node] ) 

if  abs (dualval)  >  le-08: 

nodepi [node, iteration]  +=  P [x_tilde ] *dualval 
for  (i,j)  in  vuln: 

for  defense  in  defenses: 
dualval= 

instance . dual . getValue (instance . EdgeCapacity [ i, j , defense] ) 

if  abs (dualval)  >  le-08: 

edgemu [ i, j , defense, iteration]  +=  P [x_tilde] *dualval 


expected^cost  +=  P[x_tilde]  *  obj ective_value 
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cumul_p  +=  P[x_tilde] 

#provide  intermittent  updates 
#if  scennum  %  25  ==  0: 

#  print  str(scennum)  +  ':  '  +  str (obj ective_value)  t  ' 

approx,  expected  cost:  '  t  str (expected_cost/cumul_p) 

#if  len (x_tilde) ==0 : 

#  print  '  ***'  +  str (scen_num)  +  ':  '  +  str (obj ective_value)  t 

'  approx,  expected  cost:  '  t  str (expected_cost/cumul_p) 

print  'Subproblem  Iteration:'  +  str  (iteration) 

print  '  Expected  Cost  =  '  t  str (expected_cost ) 

print  '  Worst-Case  Cost  =  '  +  str (worst_obj ) 

for  (i,j)  in  worst  seen: 

print  '  Attack:  ('  +  str(i)  +  ',  '  +  str  ( j )  +  ')  ' 

print  '  Weighted  Cost  =  '  t  str ( stoch_wt*expected_cost  t 

det_wt*worst_ob j ) 

#print  edgemu 

#print  nodepi 

if  iteration==0 : 

incumbent  =  w_hat 

UB  =  stoch_wt*expected_cost  t  det_wt*worst_obj 
inc_iteration  =  0 

print  ' ***new  incumbent:  (UB)  '  +  str (UB) 
x  star=worst_scen 
x_star_obj=  worst_obj 
x_tilde_ob j  =expected_cost 
else : 

if  UB  >  stoch_wt*expected_cost  t  det_wt*worst_obj : 

UB  =  stoch_wt*expected_cost  +  det_wt*worst_obj 

incumbent=w_hat 

inc  iteration  =  iteration 

print  ' ***new  incumbent:  (UB)  '  +  str (UB) 
x_star=worst_scen 
x_star_obj=  worst_obj 
x_tilde_obj  =expected_cost 
if  UB-LB<master_tol*LB: 
break 

for  (i,j)  in  vuln: 

if  (i,j)  in  worstscen: 

xhats [i, j , iteration]  =  1 
else : 
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xhats [i, j , iteration]  =  0 


model,  instance  =  combo_model . BuildComboMaster (Grph,  vuln 
defenses,  3,  iteration+1,  det_wt,  stoch_wt,  xhats,  nodepi,  edgemu,  supply 
capacity) 


results  =  combo_model . SolveModel ( instance ) 
instance . load (results) 

print  'Master  Iteration:  '  +  str  ( iteration) 

w_hat  =  { } 

for  (i,j)  in  vuln: 

for  d  in  defenses: 


if  instance .W [i, j , d] . value  >  le-08: 
w_hat [i, j , d]  =  1 
if  d== ' 1 ' : 

print  'Defend  ('  +  str(i)  +  ',  '  +  str  ( j )  +  ')  W= 

+  str (w_hat [i, j , d] ) 


else : 

w_hat[i,j,d]  =  0 

if  LB  <  instance . MasterObj ective () : 

LB  =  instance .MasterObj ective  ( ) 

print  'Master  objective:  (LB)  '  +  str ( instance . MasterObj ective () ) 
if  UB-LB<master_tol*LB: 
break 


print 

if  UB-LB<master_tol*LB: 

print  'Optimality  Gap  Reached  at  iteration='  +str  (iteration) 
else : 


print  'Iteration  limit  reached  at  iteration='  +str ( iteration) 
print  'Incumbent  Defense  found  at  iteration='  +str ( inc_iteration) 
for  (i,j)  in  vuln: 

if  incumbent [ i, j ,' 1 ' ]  >  le-08: 

print  '  Defend  ( '  +  str (i)  +  ' , '  +  str ( j )  +  ' ) ' 
f ile . write ('(' +str (i)  +  ','  +  str(j)  +  ')'  ) 
print  'Expected  Cost:  '  +  str(UB) 

print  'Lower  Bound:  '  +  str (LB) 

if  LB>0 : 

print  'Optimality  gap:  '  +  str ( (UB-LB) /LB) 

f ile . write ( str (weight )+' , '+str(UB)+', ' +str (xstar ) + ' , ' +str (x_star_obj ) + ' , '+str( 
x_tilde_obj ) + ' \n ' ) 

file . close ( ) 

#End  combo_main . py 
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#Begin  netx_graph.py 

#Build  the  graph  using  networkx 

import  networkx  as  nx 

def  twolines (a, b, grph,  vulnerable_edges ,  double_nodes) : 

#create  the  equivalent  of  two  parallel  edges  between  nodes  a  and  b  by 
adding 

#four  artificial  nodes  and  twelve  directed  arcs 
suf fix_l  =  "  "  +  a  +  "  "  +  b 
suf fix_2  =  "  "  +  b  +  +  a 

grph . add_node (a  +  suffix_l,  node_size  =  99, node_color=' y ' ) 

grph . add_node (b  +  suffix_l,  color='y') 

grph . add_node (a  +  suffix_2,  color='y') 

grph . add_node (b  +  suffix_2,  color='y') 

grph. add_edge (a,  a  +  suffix_l,  cost  =0) 

grph. add_edge (a,  a  +  suffix_2,  cost  =0) 

grph. add_edge (a  +  suffix_l,  a,  cost  =0) 

grph. add_edge (a  +  suffix_2,  a,  cost  =0) 

grph. add_edge (b,  b  +  suffix_l,  cost=0) 

grph. add_edge (b,  b  +  suffix_2,  cost=0) 

grph. add_edge (b  +  suffix_l,  b,  cost=0) 

grph. add_edge (b  +  suffix_2,  b,  cost=0) 

grph . add_edge (a  +  suffix_l,  b  +  suffix_l,  cost  =  1) 

grph . add_edge (a  +  suffix_2,  b  +  suffix_2,  cost  =  1) 

grph . add_edge (b  +  suffix_l,  a+suffix_l,  cost  =  1) 

grph . add_edge (b  +  suffix_2,  a+suffix_2,  cost  =  1) 

vulnerable_edges . append (  (a+suffix_l,b+suffix_l)  ) 

vulnerable_edges . append (  (a+suffix_2,b+suffix_2)  ) 

double_nodes . append (a+suf fix_l ) 

double_nodes . append (a+suf fix_2 ) 

double_nodes . append (b+suf fix_l) 

double_nodes . append (b+suf fix_2) 

return 

def  doubs  (a,  b,  grph,  vulnerable_edges)  : 

fcreate  two  directed  arcs  to  represent  an  undirected  edge  between  nodes  a 
and  b 

grph. add_edge (a,  b,  cost=l) 
grph. add_edge (b,  a,  cost=l) 
vulnerable_edges . append (  (a,b)  ) 
return 

#Create  the  model 

def  newmodel (nodes,  start_nodes,  edges,  double_edges) : 
vulnerable_edges  =  [] 
double_nodes  =  [] 
g  =  nx.DiGraph() 

#node  list 
for  i  in  nodes: 

if  i  in  start_nodes: 

g . add_node ( i,  node_color  ='b') 
else : 

g.add_node (i,  node_color='r' ) 

#edge  list 

for  i,j  in  edges: 

doubs (i, j , g, vulnerable_edges) 
for  i,j  in  double_edges : 

twolines ( i, j , g, vulnerable^edges ,  double_nodes) 
return  g,  vulnerable_edges,  double^nodes 
#End  netx_graph . py 
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#Begin  phat_code.py 
#Import  random  functions 
import  random 

from  itertools  import  combinations 
import  math 

#Take  a  list  of  vulnerable  edges  and  assign  probabilities  to  them 
def  phat_edges (graph, vulnerable^edges) : 
random. seed (0) 
phat  =  { } 
xhat  =  { } 

#Initialize  probablity  of  arcs  with  random  and  initialize  xhat  =0 
#for  tups  in  graph . edges () : 
for  tups  in  vulnerable_edges : 
temp  =  random . random ( ) 
phat [ tups ]=  temp*. 3  +  .01 
xhat  [tups]=  0 
return  phat,  xhat 

def  create_Prob (graph, vulnerable_edges, phat , max_losses=None) : 

"""Build  dictionary  of  events  and  their  associated  probabilities. 
Assumes  elements  (vulnerable_edges)  fail  independently  with 
probabilities  specified  by  the  dictionary  phat.  If 
maxlosses  is  specified,  computes  conditional  probabilities 
for  events  whose  cardinality  does  not  exceed  it.""" 

Prob={ } 

if  maxlosses  ==  None  or  maxlosses  >  len (vulnerable_edges) : 

max_losses  =  len (vulnerable_edges) 
for  i  in  range (max_losses  +  l )  : 

for  x_tilde  in  combinations (vulnerable_edges, i) : 
p=l 

for  e  in  vulnerable_edges : 
if  e  in  x_tilde: 

p  *=  phat [e] 
else : 

p  *=  l-phat[e] 

Prob [x_tilde] =p 
#normalize  probabilities 
totprob  =  math. f sum (Prob. values () ) 
for  x_tilde  in  Prob: 

Prob [x_tilde]  /  =  totprob 
return  Prob 


def  generate_xhat (vulnerable_edges, phat,  xhat): 
for  edges  in  vulnerable_edges : 
temp  =  random . random ( ) 
if  temp  <=  phat [edges]: 

xhat [edges] =1 
else : 

xhat [edges] =0 
return  xhat 

def  generate_xhat_f rom_x_tilde ( vulnerable_edges , x_tilde, xhat) : 
for  edge  in  vulnerable_edges : 
if  edge  in  x_tilde: 

xhat [edge] =1 
else : 

xhat [edge] =0 
return  xhat 
#End  phat_code.py 
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#Begin  combo_model . py 

from  _ future  import  division 

import  sys 

#Import  coopr  and  Solver 

from  coopr. pyomo  import  * 

from  coopr. opt  import  SolverFactory 

#Master  cuts 

def  MasterStochCut (model, iter_num) : 

nodeterm  =  sum (model . b [ i ] *model . pi [ i , iter_num]  for  i  in  model. Nodes) 
arcterm  =  sum (model . u [ i , j ] *model . mu [ i , j , d, iter_num] *model . W [ i , j , d]  for 
(i,j,d)  in  model .Vuln*model . Defenses) 

return  model . Z_Stoch  >=  nodeterm  +  arcterm 

def  MasterDetCut (model , iter_num) : 
edge_sum  =  ( 

sum ( (model . Y [ i, j , 'O', iter_num] * (model . c [i, j ] + (model . xhat [ i, j , iter_num] *model . q [ 
i,j])))for  (i,j)  in  model. Vuln) 

+sum ( (model . Y [ j , i , 'O', iter_num] * (model . c [ j , i ] + (model . xhat [ i , j , iter_num] *model . q 
[j,i])))for  (i,j)  in  model. Vuln) 

+sum ( (model .Y [i, j ,' 1 ', iter_num] * (model . c [i, j ]) )  for  (i,j)  in 

model .Arcs) 

) 

node_sum  =  sum (model . v [ j ] *model . S [j ,  iter_num]  for  j  in  model. Nodes) 
return  model . Z_Det  >=  edge_sum  +  node_sum 

def  MasterDetFlowBalance (model, node, iter_num) : 

amountln  =  sum  (model . Y [ i , j , def ense, iter_num]  for  ( i, j , defense)  in 
model . Arcs*model . Defenses  if  j==node) 

amountOut  =  sum  (model . Y [ i, j , defense, iter_num]  for  ( i , j , def ense)  in 
model .Arcs*model . Defenses  if  i==node) 
supply  =  model. b[node] 

return  amountOut  -  amountln  -  model . S [node, iter_num]  <=  supply 

def  MasterDetEdgeCapacity (model, i, j , defense, iter_num) : 

return  model . Y [ i, j , defense, iter_num]  +  model . Y [j , i, defense, iter_num]  <= 
model . u [i, j ] *model . W [i, j , defense] 

#Choose  one  defense  for  each  vulnerable  edge 
def  MasterOneDefense (model,  i,  j )  : 

return  sum (model . W [ i, j , d]  for  d  in  model . Defenses )  ==  1 

#Limit  number  of  defenses 
def  MasterDef enseLimit (model ) : 

return  sum (model . W [ i, j ,' 1 ' ]  for  (i,j)  in  model. Vuln)  <=  model ,max_defenses 

#Master  Objective 

def  MasterObj ective (model ) : 

return  model . det_wt*model . Z_Det  +  model . stoch_wt*model . Z_Stoch 

#Ensure  name  of  networkx  graph  is  Grph  to  run  effectively 
#  need  networkx  graph  named  Grph,  dictionary  for  xhats 
def  BuildComboMaster (Grph,  vulnerable_edges ,  defenses,  max_defenses, 
iterations,  det_wt,  stoch^wt,  xhats,  nodepi,  edgemu,  supply,  capacity) : 

#Define  model  type 
model  =  AbstractModel ( ) 

#Put  model  into  Pyomo 

nodes  =  Grph. nodes () 

model. Nodes  =  Set (initialize=nodes) 
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model. Vuln  =  Set (within=model . Nodes*model . Nodes , 
initialize=vulnerable_edges ) 
arcs  =  Grph. edges () 

model .Arcs  =  Set (within=model .Nodes*model .Nodes, initialize=arcs) 

iters  =  range (iterations) 

model. Iters  =  Set (initialize=iters) 

model . Defenses  =  Set ( initialize=defenses) 

model .max_defenses  =  Param (initialize=max_defenses) 

#Node  dual  parameters 

model. pi  =  Param (model . Nodes*model . Iters,  initialize=nodepi ) 

#Arc  capacity  dual  parameters 

model. mu  =  Param (model .Arcs*model . Def enses*model . Iters,  initialize=edgemu) 

model. xhat  =  Param (model .Vuln*model . Iters,  initialize=xhats) 

♦arc  costs 
values=  { } 

for  item  in  Grph . edges () : 

values[item]  =  Grph . edge [ item [0 ]] [item [ 1 ]][' cost ' ] 
model . c  =  Param (model .Nodes, model .Nodes,  initialize=values) 

#supplies 

model. b  =  Param (model .Nodes, initialize=supply) 

♦capacities 

model. u  =  Param (model .Vuln,  initialize=capacity) 

♦Penalty  model 

model. q  =  Param (model .Arcs,  initialize=100) 

♦Penalty  unsupplied  nodes 

model. v  =  Param (model .Nodes,  initialize=10) 

model. det_wt  =  Param ( initialize=det_wt) 
model . stoch_wt  =  Param (initialize=stoch_wt) 

♦define  objective  variables 

model . Z_Stoch  =  Var (domain=NonNegativeReals ) 
model . Z_Det  =  Var (domain=NonNegativeReals ) 

♦define  design  variables 

model. W  =  Var (model .Vuln*model . Defenses,  domain=Binary) 

♦define  response  variables  for  deterministic  cuts 
model. Y  =  Var (model .Arcs*model . Defenses*model . Iters, 
domain=NonNegativeReals ) 

♦Define  fuel  shortfall  at  Node 

model . S  =  Var (model .Nodes*model . Iters,  domain=NonNegativeReals) 

♦Master  cut  constraints 

model .MasterStochCut  =  Constraint (model . Iters,  rule=MasterStochCut) 
model .MasterDetCut  =  Constraint (model . Iters,  rule=MasterDetCut) 

♦Det  flow  constraints  for  each  iter 

model .MasterDetFlowBalance  =  Constraint (model .Nodes*model . Iters, 
rule=MasterDetFlowBalance) 
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model .MasterDetEdgeCapacity  = 

Constraint (model . Vuln*model . Defense s*model . Iters,  rule=MasterDetEdgeCapacity) 
#Master  arc  defense  constraints 

model .MasterOneDefense  =  Constraint (model .Vuln,  rule=MasterOneDefense) 
#Master  limit  on  number  of  defenses 

model .MasterDefenseLimit  =  Constraint (rule=MasterDefenseLimit) 

#Objective  function 

model .MasterObj ective  =  Objective (rule=MasterObjective,  sense=minimize) 
instance=model . create  ( ) 
return  model,  instance 

#Create  the  model  instance  and  Solve  returning  values  for 

#min  flow  path  in  variable  called  shortest 

#need  list  of  double_nodes  to  create  key  and  edge  labels 

def  SolveModel ( instance) : 

opt  =  SolverFactory ( "cplex" ) 

#opt  =  SolverFactory ( "glpk" ) 

#ask  for  duals  and  reduced  costs  from  solver 
results  =  opt . solve (instance,  suf fixes= [ ' dual ' , ' rc ' ] ) 
return  results 

def  BuildResults (instance, results) : 
instance. load (results) 

#figure  arcs  defended 

defended  =  [] 

for  key  in  instance. W: 

if  instance . W [ key ]. value  >  0: 
defended. append (key) 

obj ective_value  =  instance. MasterObjective ( ) 
return  defended,  ob j ective_value 
#End  combo_model . py 
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#Begin  pyo_model.py 

from  _ future  import  division 

import  sys 

#Import  coopr  and  Solver 

from  coopr. pyomo  import  * 

from  coopr. opt  import  SolverFactory 

#Network  balance  of  flow  constraint 
def  FlowBalance (model , node)  : 

amountln  =  sum  (model . Y [i , j , def ense ]  for  (i, j , defense)  in 
model .Arcs*model . Defenses  if  j==node) 

amountOut  =  sum  (model . Y [ i, j , defense]  for  ( i, j , defense)  in 
model .Arcs*model . Defenses  if  i==node) 
supply  =  model. b[node] 

return  amountOut  -  amountln  -  model . S [node]  <=  supply 

#Define  objective  equation 
def  FlowCost (model) : 
edge_sum  =  ( 

sum ( (model . Y [i, j ,  '0']* (model . c [ i, j ]  +  (model . xhat [i , j ] *model . q [ i, j ] ) ) ) for  (i,  j ) 
in  model. Vuln) 

+sum ( (model . Y [ j , i, '0']* (model . c [ j , i ] + (model . xhat [ i, j ] *model . q [ j , i ] ) ) ) for  (i, j ) 
in  model. Vuln) 

+sum ( (model . Y [i, j 1 ']* (model . c [i, j ]) )  for  (i,j)  in  model. Arcs) 

) 

node_sum  =  sum (model . v [ j ] *model . S [j ]  for  j  in  model. Nodes) 
return  edge_sum  +  node_sum 

def  EdgeCapacity (model,  i,  j,  defense): 

return  model . Y [i, j , defense]  +  model . Y [j , i , def ense ]  <= 
model . u [i , j ] *model . w_hat [ i, j , defense] 

#Ensure  name  of  networkx  graph  is  Grph  to  run  effectively 
#  need  networkx  graph  named  Grph,  dictionary  for  xhats 
def  build_model (Grph,  xhat,  w_hat,  vulnerable_edges ,  defenses,  supply, 
capacity) : 

#Define  model  type 
model  =  AbstractModel  ( ) 

#Put  model  into  Pyomo 
nodes  =  Grph. nodes () 

model. Nodes  =  Set (initialize=  nodes) 
edges  =  Grph. edges () 

model. Arcs  =  Set (within=model .Nodes*model .Nodes, initialize  =  edges) 

model. Vuln  =  Set (within=model .Nodes*model .Nodes, initialize  = 
vulnerable^edges) 

model . Def enses  =  Set ( initialize=defenses) 

fcreate  suffixes  for  dual  information 

model. dual  =  Suffix (direction=Suf fix. IMPORT_EXPORT) 

#Create  model  for  costs 

values=  { } 

for  item  in  edges: 

values[item]  =  Grph . edge [ item [0 ]] [item [ 1 ]][' cost ' ] 
model . c  =  Param (model .Nodes, model .Nodes,  initialize=values) 
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#Supplies 

model. b  =  Param (model .Nodes, initialize=supply) 

#Penalty  model 

model. q  =  Param (model .Arcs,  initialize=100) 

#Penalty  unsupplied  nodes 

model. v  =  Param (model .Nodes,  initialize=10) 

#Capacities  on  vulnerable  arcs 

model. u  =  Param (model .Vuln,  initialize=capacity) 

#Define  attack  variable 

model. xhat  =  Param (model .Vuln, initialize=xhat) 

#Define  attack  variable 

model. w_hat  =  Param (model .Vuln*model . Defenses, initialize=w_hat) 
#Define  Variable  we  want  to  change 

model. Y  =  Var (model .Arcs*model . Defenses,  domain=NonNegativeReals) 
#Define  fuel  shortfall  at  Node 

model . S  =  Var (model .Nodes,  domain=NonNegativeReals) 

#Constraint  for  balance 

model . FlowBalance  =  Constraint (model . Nodes,  rule=FlowBalance) 

model . EdgeCapacity  =  Constraint (model .Vuln*model . Defenses, 
rule=EdgeCapacity) 

#Objective  function 

model . FlowCost  =  Objective (rule=FlowCost,  sense=minimize) 
instance=model . create ( ) 

return  model,  instance 

#Create  the  model  instance  and  Solve  returning  values  for 

#min  flow  path  in  variable  called  shortest 

#need  list  of  double^nodes  to  create  key  and  edge  labels 

def  solve_model (instance) : 

opt  =  SolverFactory ( "cplex" ) 

#opt  =  SolverFactory ( "glpk" ) 

#ask  for  duals  and  reduced  costs  from  solver 
results  =  opt . solve (instance,  suf fixes= [ ' dual ' ,  ' rc ' ] ) 
return  results 

def  build_results (instance, results)  : 
instance. load (results) 

#figure  arcs  used 
shortest  =  [] 
edge_values  =  { } 
for  key  in  instance. Y: 

if  instance . Y [ key ]. value  >  0: 

shortest . append ( ( key [ 0] , key [1 ] ) ) 

edge^values [ (key[0] , key [1] ) ]  =  instance . Y [key ]. value 

ob j ective_value  =  instance . FlowCost (instance) 
return  shortest,  edge_values,  ob j ective_value 

#End  pyo_model.py 
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